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ABSTRACT 

In this paper we will consider the problem of the numerical simulation of 
non-Gaussian, scalar random fields with a prescribed correlation structure pro- 
vided either by a theoretical model or computed on a set of observational data. 
Although, the numerical generation of a generic, non-Gaussian random field is 
a trivial operation, the task becomes tough when constraining the field with a 
prefixed correlation structure. At this regards, three numerical methods, useful 
for astronomical applications, are presented. The limits and capabilities of each 
method are discussed and the pseudo-codes describing the numerical implemen- 
tation are provided for two of them. 

Subject headings: methods: data analysis - methods: numerical 

1. INTRODUCTION 

Computer-aided modeling is becoming an essential tool in designing new experiments 
and in testing theoretical models against the observational data. For example, because of 
the cost of any space-based telescope, nowadays it is not even conceivable to plan a mission 
without first simulating the performances either of the instruments and/or of the observing 
mode. It is obvious to stress that the reliability of such simulations depends critically on 
the possibility to reproduce realistic physical scenarios. 

A wide assumption, which is often made because of its simplicity, is that the processes 
underlying a given physical phenomenon obey Gaussian, and therefore linear, statistics. 
However, although as practical as this assumption could be, it is not applicable for most of 
the physical systems which, on the contrary, are expected to be characterized by nonlinear 
behaviours. Some examples: 
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- the high spatial resolution observations of sky images are revealing a lot of details of 
the sky emission which are not of easy interpretation (e.g. Herbstmeier et al. (1998)). 
For istance, in the far infrared spectral domain, the studies of source properties imply 
the disentangling of the source emission and position from those of a much stronger 
background whose spatial structure is highly non-Gaussian. This task becomes crucial 
for most space-borne surveys of the extragalactic sky and of star-forming regions in 
the Galaxy for which it is of interest to simulate the far-IR and sub-mm emission 
of the Galaxy and the source confusion in the beam. Their nature is intrinsically 
non-Gaussian and to match their observed properties an appropriate method must be 
used. 

- The interpretation of the new flow of data on the Cosmic Microwave Background 
(CMB) spatial distribution, from present and future experiments (BOOMERANG, 
MAXIMA, MAP and PLANCK), is involving a large theoretical effort (see e.g. Verde 
et al. (2000); Matarrese, Verde & Jimenez (2000); Contaldi & Magueijo (2001) and 
references therein). Studies of the spatial structure of the CMB provide fundamental 
clue on the physical processes generating the primordial density fluctuations which are 
thought to be at the origin of the present-day structures. There are deep theoretical 
motivations, both in the framework of inflationary models and in cosmological 
defects scenarios, to consider the initial density perturbations obeying non-Gaussian 
statistics. In any case, the subsequent growth of these density fluctuations, triggered 
by the gravitational potential, makes the late time evolution nonlinear. The testing of 
these predictions requires algorithms able to discerne the true statistical nature of the 
observed fields. Some of their properties can be analytically recovered but the bulk of 
non-Gaussian random fields characteristics can be inferred only through simulations 
(see e.g., Moscardini et al. (1991)). 
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The aim of this paper is to provide a general and mathematical approach to the problem 
of generating non-Gaussian random fields when a correlation structure is given either by 
a theoretical model or by the statistical analysis of experimental data. Some numerical 
procedures to perform simulations of such fields are also presented. The arguments are 
outlined on a quite general basis in view of applications in a wide astrophysical context. The 
reason to fix the correlation structure is that this is the simplest way to obtain non-trivial 
(i.e. non-pure noise) fields (see below). 

The problem can be formalized as follows. A real, random field R(t) can be defined 
as a collection of random variables {r} at points with coordinates {t} = {(ti,t 2 , ■ ■ ■ ,t n )} 3 
in a n-th dimensional "parameter space". In other words, for each "position" t, R(t) = r, 
where r is a vector characterized by a multidimensional distribution function Fr(t) and a 
multidimensional probability density function /i?(r). According to the particular problem 
at hand, {t} may correspond to a set of spatial/angular coordinates (spatial random fields), 
to time (time processes), to a mix of these two (spatio-temporal random fields) or even to 
more general situations. 

In many practical situations they are of interest the so called scalar random fields, 
where r = r. This means that for a specific t, R(t) is characterized by a scalar random 
variable r with one-dimensional distribution function (DF) F R (r) and one-dimensional 
probability density function (PDF) fnij-). In this paper we will consider only this kind of 
random fields. The more general case of vector-valued random fields represents a more 
complex problem and will not be addressed in this work (for more details about this topic 
see Popescu, Deodatis & Prevost (1998)). 

3 From now on, in order to distinguish them from scalar quantities, we will denote vector 
quantities in boldface. 
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The main problem in the numerical simulation of a generic R(t) is that, in general, 
given two arbitrary "positions", say ti and t 2 , R(ti) and R(t 2 ) are not independent 
one of the other. As well known from the standard theory of random processes (e.g. 
Grigoriu (1995)), the practical applications to generate random fields requires to put some 
constraints. The most common choice is that, specified the distribution function Fn(r), 
R(t) be completely characterized by the covariance function, ^ R (ti,t 2 ), with 



where E[.] stands for expected value. The reason is that ^(ti,t 2 ) represents the simplest 
form of mutual relationship between the elements of R(t). 

In astronomical applications, very often it is possible to adopt some simplifying 
conditions. In particular, it is possible to assume that R(t) is isotropic. This means that 
the covariance function depends on the length of the vector t x — t 2 but not on its direction: 
{i?(ti,t 2 ) = G?(l|ti — t 2 ||) 4 . In other words, R(t) is characterized by a spherical symmetry. 
This property is very useful since it allows to characterize R(t) through the correlation 
function 



where r = ||r||, and and a\ are, respectively, the mean and the variance corresponding 
to the distribution F R (r). 

Although the isotropic case is of large interest in astronomical applications, here we 
prefer to adopt a more general formalism, suited for all homogeneous fields. Indeed, in 
this case, the covariance function depends on t x — t 2 . According to this definition, p(r), in 
equation (2) has to be replaced by p(r). 

4 We remind that for a column vector a = (ai, a 2 , . . . , a7v) T , ||a|| = [a 7 ^] 1 / 2 = (J2iLi a ?) 1//2 
provides its length (norm). Here a T means the transpose of vector a. 



^(t 1 ,t 2 ) = E[i?(t 1 ) i?(t 2 )], 



(1) 



Pr(t) = E 



(R(t) - fj, R ) (Rjt + t) - /ir) 



(2) 
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2. 



PRELIMINARY NOTES 



Most of the techniques for simulating a non-Gaussian, scalar, random field R(t), with 
a prescribed correlation function, Pr(t), and a prescribed one-dimensional marginal F R (r), 
are more or less explicitly based on the two following steps: 

- generation of a zero-mean, unit- variance, scalar, Gaussian random field X(t) with a 
prefixed correlation structure pxij); 

- mapping (transformation) X(t) — > R(t) according to 



where g[.} represents an appropriate function. This operation is named memoryless 



The rationale behind such a procedure is that the direct generation of a generic R(t), 
with a specific p r (r), is a very difficult operation. The Gaussian case represents a useful 
exception. Hence, it results much easier to obtain R(t) by transforming a precomputed 
X(t). However, after the mapping (3), in general Px(t) does not coincide with Pr{t). 
Therefore, it is necessary to transform an X(t) characterized by an appropriate pxir) 
whose functional form depends on Pr{t). 

It is well known from elementary statistics that for a mapping r = g(x), with x the 
standard one-dimensional Gaussian variable, the PDF of the random variable r can be 
obtained from that of the random variable x via a change of variable technique. In the 
general case that the tranformation g(.) is not one-to-one, if the equation 



R(t) = g[X(t)l 



(3) 





(4) 
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has a numerable set of M real solutions {xi(r),X2(r), . . . , xm(j')}, and if g'j = [dg(x)/dx] x=Xj , 
j = 1, 2, . . . , M exist, then fn(r) is given by (Papoulis 1991) 

In correspondence to the values r*, where equation (4) does not have real solutions, it 
happens that /jj(r*) = 0. Furthermore, the correlation function Pr(t) is given by (Grigoriu 
1995) 

Pr{t) = — / [ff(xi) - [#(z 2 ) - 0(xi,X2;px(t)) ttei cfc 2 , (6) 

® R J -oo J -oo 

where, x\ = x(t) and x 2 = x(t + r), and 

1 / arf + a;| -2p x (r):ri:r 2 \ 
^^(r)) = 27r(1 _ p|(r))1/2 exp ^ 2(1 _^ (T)) j • ( 7 ) 

At first sight, from these equations it may seem that, given the appropriate function 
g(.) and the covariance function px(t), obtained via the inversion of equation (6), it is 
possible to generate an R(t) characterized by an arbitrary Pr(t). In reality, given a generic 
(?(.), there is no guarantee that equation (6) can be inverted. Furthermore, it is possible to 
show (Ogorodnikov & Prigarin 1996) that Pr{t) can take values only in the interval 

Pi? (T)e[ P *,l], (8) 

where 

P* = (jf Fr 1 ^) - «) d« - /^) , (9) 

with 

F R \a) = inf{r : F R (r) > a} (10) 

providing the smallest value of the random variable r satisfying the condition that 
Fr( t ) > a - In particular, p* = —1 only for symmetric distributions. For example, 
p* ~ —0.645 in case of the exponential PDF: /ij(r) = /3 exp(— (3r),r > 0. This shows that, 
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in general, for a fixed g(.) it is not possible to obtain a p* R {r) presenting values external to 
the interval (8). 

Another problem stems from the fact that the function px(t), necessary to obtain the 
target Pr{t) via equation (6), must be a non-negative definite function 5 otherwise the 
generation of the aimed R(t) is prevented since px(j) does not represent any correlation 
function. 

In principle, equation (6) can be used to obtain Px(t) in a closed form, but in reality 
such an approach can be followed only in a limited number of cases. Indeed, very often 
the transformation (6) has a very complex form and can be inverted only via a numerical 
approach. 

3. ANALYTICAL METHOD 

Certainly one of the most effective method for generating R(t) is represented the 
analytical handling of equations (5) and (6). Unfortunately, this is also the most difficult 
approach to pursue; only in a limited number of cases it has been possible to find out the 
analytical relationship between pxir) and Pr{t). Some useful examples are presented 
below (see also figure 1): 

5 It should be remembered that only for a non-negative defined function the corresponding 
Fourier transform has non-negative values. Therefore, pxij) must share this property since, 
according to the Wiener-Khinchin theorem, the power-spectrum of a process is provided by 
the Fourier transform of the corresponding correlation function. 
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3.1. Lognormal Fields 

If X(t) is a homogeneous, zero-mean, unit- variance Gaussian fields with a correlation 
function Px(t), the random fields obtained via 

L(t) = e» +aX V (11) 

are called Lognormal fields since they are characterized by the one-dimensional marginal 
Lognormal PDF 

f L (l) = 1 e -(lni-/,) 2 /(2 CT 2 ) ; I > 0. (12) 

I ay/271 

It is possible to show (Vanmarke 1984) that the moments of order k of L(t) are given by 

E[L k ] = e k » +k2 ° 2 / 2 . (13) 

In particular, the mean and the variance are 

= e^ 2 / 2 , a\ = e^ +a2 (e^ - 1), (14) 

respectively Furthermore, it is possible to show that the relationship between Pl(t) and 
Px(t) can be expressed in the form 

a 2 px(T) _ , 

Pl(t) = - 2 - ■ (15) 
e — 1 

From this equation it is trivial to see that, when a = 1, the lower bound of Pl(t) is 
~ -0.368. 



3.2. Gamma Fields 

If X s (t), s = 1,2, ... ,2m, is a collection of independent zero-mean, unit-variance 
Gaussian fields with the same correlation function p x (t), the random fields obtained via 

2m 

G m (t) = -$> s 2 (t) (16) 

s=l 
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are called Gamma fields. That is because the corresponding one-dimensional marginal PDF 
is a Gamma distribution with m degrees of freedom 

fcjg) = f ^- ) g m - 1 e- 9 , g>0, (17) 

where T(.) is the Gamma function. 

It can be shown that the moments of order k of G m are given by 

E[G ™ ]=£ r^ = n (m+2) ' k> - m - (i8) 

^ ' i=0 
In particular, the mean and the variance are 

Pc m = cr Gm = m. (19) 

It can also be shown (Hasofer, Ditlevsen & Tarp-Johansen 1998) that, idependently from 
the value of m, the relationship between PG m ( T ) an d Px(t) can be expressed in the form 

PgJt)= P 2 x (t). (20) 

The class of the Gamma fields is interesting since it contains, as particular cases, both the 
Chi-Square and the Exponential fields. 

From equation (20) it appears that the lower bound of PG m ( T ) is zero - m other words, 
through the mapping (16) it is not possible to obtain G m (t) characterized by correlation 
functions with negative values. Here, however, it is necessary to stress that such a limit is 
not intrinsic to the Gamma fields, but only to the transformation (16). In fact, through 
different mappings it is possible to generate G m (t) with correlation functions having 
negative values (see below). 
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3.3. Beta Fields 

Given two independent Gamma fields, say G m (t) and G n (t), characterized by the same 
correlation function Pg(t), the random fields obtained via 

B - {t) = cJnkit) (21) 

are called Beta fields because their one-dimensional marginal PDF is a Beta(m, n) 
distribution 

f Bm M = j^—, x m -\l - x) n ~\ < b < 1. (22) 
B[m, n) 

It can be shown that the moments of order k of B mn are given by 

k _ T(m + k)T(m + n) 
HB l ~ T(m)T(m + n + ky 1 ' 

In particular, the mean and the variance are 

m 2 mn 

m + n (m + n) z (m + n + 1) 

respectively. 

It can be also shown (Hasofer, Ditlevsen & Tarp-Johansen 1998) that the relationship 
between PB mn (r) and px(r) can be expressed in the form 

pB mn (r ) = 1 - S m+n [p x {r)}, n + m>l, (25) 

where 

/ 1 — n \ q r 1 / — n \ 

0< P <1, gG{l,2,...}, (26) 

with the end values 5^(0) = 1 and S^l) = 0. The class of the Beta fields is basic for 
describing variables bounded at both sides. For example, -Bn(t) corresponds to the Uniform 
field. 

As well as for the Gamma fields, also for the Beta fields obtained through mapping 
(21) it happens that the lower bound of pB mn (. T ) lii zero. Again, this limit is not intrinsic to 
the Beta fields but only to the particular mapping used. 



w .,(i^)*U_ rt _£i(^) 
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4. 



NUMERICAL METHOD 



In case one is interested in a R(t) characterized by an Fn(r) not reproducible through 
a function g(.) available in analytical form and/or that does not permit an easy calculation 
of px{t), it is necessary to resort to numerical methods. Regarding this, the following 
presents three methods that can be useful in astronomical applications. 



The most obvious method is based on the numerical inversion of equation (6). 
Whenever possible, this is the 'method' to use since, contrary to the procedures presented 
below, it is able to provide exact results (within the limits of the numerical computation). 
In particular, this inversion operation is feasible when g(.) is a monotonic increasing, real 
function. Indeed, in this case the relationship between Pr{t) and px(r) can always be 
inverted. 

In the case of distributions F r with no atoms (a concentration of a finite probability 
mass at a point), a very useful kind of monotonic increasing functions g(.) is represented by 
the mapping 



where Fx denotes the Gaussian distribution function and F R the inverse distribution 
function of R(t). Indeed, through g(.) = F R 1 {F X (.)} the generation of fields -R(t), with 
arbitrary one-dimensional marginal distribution functions, is possible. Furthermore, it can 
be shown (Ogorodnikov & Prigarin 1996; Grigoriu 1995) that via the mapping (27) it is 
possible to obtain R(t) characterized by Pr{t) fully exploiting the interval (8) with p* that 
can be simplified to the form 



4.1. Change of Variable Method 



R(t) = F R 1 {F x [X(t)}}, 



(27) 



E[g{x) g{-x)}-p 2 R 
4 



(28) 



Figure 2 shows the relationship between px and the pr, concerning some well known 
F R , obtained via equations (27) and (6). From this figure it is possible to realize some 
interesting points that can also be proved via theoretical arguments (Grigoriu 1995, 1998): 



- Pr(t) is an increasing function of pxij); 



Pr(t)\ < \ Px (t)\; 



- the difference between Pr(t) and pxij) are not significant for a broad range of values 
of these functions. 

As explained in Section 2, once Px{t) has been calculated, it is necessary to check that this 
function is non-negative definite. A possibility consists in evaluating the Fourier transform 
of Px{t) and verifying that it presents no negative values. 

The only concern regarding the numerical inversion of equation (6) is that, in general, 
this operation requires the calculation of a large number of double integrals. In certain 
situations, that could represent a computationally too expensive problem, and therefore it 
is necessary to resort to other numerical techniques. 



An alternative approach for the simulation of a non-Gaussian R(t) is based on the 
expansion of the field in Hermite polynomials 6 . These polynomials can be defined through 
the Rodriguez's formula 



6 This expansion, known also as Edgeworth expansion, was already applied in a Cosmo- 
logical context (see Colombi (1994)) 



4.2. Hermite Expansion Method 




d" 



(29) 



n = 0,l,..., 



and have the important property of being orthogonal relative to the standard Gaussian 
distribution, so that 



/+oo 1 
R m (x)R n (x) -= e- x2 ' 2 dx = n\ 5 mn (30) 
■oo V^TT 



where 5 mn is the Kronecker function. 



An explicit expression for H n (x) is given by Blinnikov & Moessner (1998) 

(-l) k x n ~ 2k 
k\{n-2k)\ 2 fc ' 



fc=0 

where |_^J means the largest integer k < z. 
A field R{t) can be expanded according to 

iT(t)= J> fe H fe pT(t)), (32) 

fc=0 

where the coefficients {a^} are unknown and must be determined. In the practical 
applications, also the "optimal" value N H has to be determined. 

One possibility for obtaining the coefficients {a^}, for a fixed N H , is to minimize the 
objective function (Grigoriu 1995) 

5 2 = E[\R(t) -iT(t)| 2 ] (33) 

that yields the conditions 



E 

so that 



N H 

R(t)-J2ai H,(X(t)) ] H fc (X(t)) 

1=0 



0, k = 0,l,...N H , (34) 



a k = ^ E[R(t) H fc (X(t))], k = 0, 1, . . . , N H , (35) 



because of equation (30). Here, the important point is that the coefficients {a^} are 
independent from the structure of X(t) since, for a specific position t, the value of R(t) 
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depends only on X(i). That allows us to estimate such coefficients by means of 

a k = i- E[r R k (x)}, (36) 
k\ 

where x is the standard Gaussian random variable, and r is the random variable distribuited 
according to the marginal distribution required for R{t). Following this approach, the 
procedure implemented in the subroutine HermCoeff in figure 3 is: 

1. generation of a large (column) array of independent and uniform random deviates 
u = [ui,u 2 , ■ ■ .,u N } T ; 

2. mapping of u in two arrays x = [xi,x 2 , ■ ■ ■ , Xn] T = F^in) and r = [r\, r 2 , . . . , r^] T = 
F i ^" 1 (u). Here, x is an array of independent standard Gaussian random deviates, 
whereas r is an array of independent random deviates distributed according to F R (r); 

3. calculation of the arrays h k = h k 2, ■ ■ ■ , hkN] T with h ki = H fc (xj), k — 0,1, ... , N H ; 

4. calculation of the coefficients {a k } according to 

ak = ITn rT hk - (37) 

The "optimal" value for Nh can be determined on the basis of the value of the parameter e 
provided by the criterion 

e = DIST[/ Jlj / r .] J (38) 

where DIST[., .] is a measure of the distance between }'r{j-) and the PDF, f r *, relative to 
the random deviates 

N H 

r* = ^a fc H fc (x). (39) 

k=0 

Although this approach also presents the problem that Pr{t) ^ Px(t), here the 
situation is easier than the method considered in the previous section. Indeed, because of 
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the orthogonality of the Hermite polynomials (Grad 1949), and in particular because of the 
so called Kibble-Slepian formula (Slepian 1972; Declercq 1998), we have that (Declercq 
1998; Sakamoto k Ghanem 1999) 

PR{T) Nh • (40) 

E fc =i kl 4 

Therefore, Px(t) can be obtained by the numerical inversion of a polynomial function. It 
is better to recall that, before using it, such px{r) must be tested to be a non- negative 
definite function. 

Once {afc}, Nh, and pxij) have been determined, R(t) can be obtained by equation 
(32) that is implemented in the subroutine Field_Herm shown in figure 3. 

Some notes on the use of the subroutine HermCoeff: 

- the input parameters are the length of the arrays x and r, the PDF /^(r), the target 
correlation function Pr(t), and the parameter e for the convergence criterion. The 
output quantities are the number Nh of terms for the Hermite expansion, the vector 
a = {g^} containing the values of the coefficients of the expansion, and the correlation 
function Px(t) of the Gaussian random field X(t) to be used in the subroutine 
Field_Herm; 

- typical value for iV is several thousands; 

- for the stopping criterion, DIST[., .] < e, it is necessary to choose a distance measure 
between /#(r) and the corresponding approximation f r *. Such a choice, as well as the 
value of the parameter e, is very situation dependent. One possibility is to calculate 
the difference between the corresponding histograms. Take note, however, that this 
method can be troublesome in case of very skewed distributions. In this case it is 
advisable to resort to the methods of PDF estimation that do not make use of binning 



-17- 



of the data as, for example, the kernel and the Johnson empirical distributions 
methods (Vio et al. 1994). 

Some conveniences concerning the algorithm: 

- in case of isotropic random fields, it is possible to work with an one-dimensional 
correlation function Pr(t); 

- once the coefficients a of the expansion have been determined, these can be used for 
simulating an unlimited number of random fields R(t); 

- the algorithm works also in case of very skewed distributions (see figure 4). 

One inconvenience in using this algorithm is that /ij(r) is only approximated and in 
particular situations this fact can be troublesome. For example, in case of strictly positive 
random fields R(t), it could happen that R*(t) presents some negative values. However, 
if the approximation is good enough (e.g. only a few values violate the constraints), the 
solution to this kind of problem can be very simple (es. the reflection of the negative values 
to positive values). 

4.3. Method of Yamazaki & Shinozuka 

Figure 5 shows the subroutine Field-ID F implementing an algorithm based on an 
idea by Yamazaki & Shinozuka (1988). The rationale behind this code is simple and is 
based on an iterative procedure. As explained in section 2, the mapping (27) deforms 
Px(t) — > Pr(t), and consequentely the corresponding power spectrum 5x(k) — > S#(k), 
in a complex way. However, if one applies the transformation (27) to an initial jW(t), 
characterized by a power spectrum S^(k) set equal to the target Sn(k), it is possible to 
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recover information on the relationship between S^(k) and Sn(k) from the power spectrum 
S^(k) of i?«(t). A new Gaussian field X( 2 )(t), characterized by a power spectrum S^(k), 
is then built with the aim that, after mapping (27), S^(k) is closer to Sn(k) than S'^(k) . 
This operation is carried out at line 21 of the code where the power spectrum of X^ +1 ^(t) 
is assumed to be given by 

4 +1) (k) = c(k)4 ) (k), (4i) 

with 

C(k) = 4t^- (42) 

Through this step, (k) is modified according to the fractional difference, C(k), between 
S#(k) and S^(k). The entire procedure can be repeated n times until that <S^(k) is a 
good approximation of Si?(k) . 

The code presented in figure 5 has been modified with respect to the original version 
of Yamazaki & Shinozuka. The main difference referes to the implementation of steps 15, 
19-20, 23-25. The task of these steps is to constrain the range of the permitted values for 
C(k). Indeed, strictly speaking, the use of such a factor is correct only within the hypothesis 
that the map (27) is linear. The consequence is that in many situations C(k) appears as a 
highly oscillating function with large extremes even in cases where the target Sr(\s), and 
therefore the starting S^(k), is a smooth function. Because of this fact, in general the 
algorithm of Yamazaki & Shinozuka converges only with moderately non-Gaussian fields. 
Our modification is based on the idea that, although some values of C(k) could be too 
large, they are still able to provide indications concerning the direction of the corrections to 
make in Sj^(k), via equation (41), for improving the results of the iterative process. That 
suggests the following procedure: 



1. once C(k) was computed, the subset k* of the indices k has to be identified for which 
C(k*) is larger than a threshold 1 + 5, where 5 is an appropriate value; 
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2. set C(k*) = 1 + 5. In this way, it is possible to obtain a smoothed version of C(k) 
that maintains the original information on the direction of the correction for each 
frequency k; 

3. if after this operation it happens that DIST[^ } (k), S R (k)\ > DIST[^ _1) (k), S R (k)}, 
where DIST[., .] indicates a distance measure between the two arguments, it is 
necessary to rescale the parameter S according to a prescribed schedule. This point 
makes it possible to avoid troublesome oscillations of the algorithm. 

From our simulations, it appears that after these modifications the algorithm converges also 
in situations where the original method fails. 

Some notes on the use of the subroutine FieldJDF: 

- the first two input quantities are the phase angles, 0(k), and the power spectrum, 
S'x(k), of a zero- mean, unit- variance, and Gaussian random field. Here, the important 
point is that Sx(t) is set equal to the target S#(k). 

The third input quantity is the initial value, Sq, of the parameter 5. Typically, 
5 = 1-2, but such a choice is not so critical for the final results. 

For the fourth input quantity, e, see below; 

- in the stopping criterion, DIST[., .] < e, any measure of distance can be used between 
<S^(k) and Sfl(k). An interesting suggestion comes from Popescu, Deodatis & Prevost 
(1998) that in their work use the quantity 

DIST[S<?(k), S R (k)] = Ek 'f ( ^r^ (k)l - (43) 
Typical value for e are of order of 10~ 2 -10" 3 ; 

- the scaling, SCALE[.], of the parameter S can follow any schedule. In our simulations 
we have halved the value whenever required by the convergence check. 
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In the context of the astronomical applications, some limitations concerning the 
algorithm are: 

- the target power spectrum Sfl(k) must be a smooth function. That means to work 
with the expected power spectra of the fields (NB. in certain engineering applications 
this is a demanded point). Such requirement is due to the fact that the generation 
of R(t), according to the procedure implemented in the algorithm of figure 5, in 
practice constitutes an optimization problem. Since the rougher a function, the larger 
is the corresponding number of degrees of freedom that must be accounted for by an 
optimization procedure, a non-smooth Sft(k) will be hardly a solvable problem with 
the present algorithm; 

- although the algorithm is more robust than the original version of Yamazaki & 
Shinozuka, it still presents convergence difficulties in case of PDFs very different from 
the Gaussian one (see figure 6). In particular, the most serious problems concerns 
very skewed distributions. The reason can be understood from the figures 7a-d, where 
the mapping (27) is presented for four distributions Xd with d = 1-4. The first two 
distributions represent situations that the algorithm is not able to solve, the third 
one corresponds to a difficult case, whereas the last distribution can be easily handled 
with. It is easy to see that the most problematic situations concern the mappings 
where a large portion of the domain of the Gaussian random variable is projected onto 
an almost constant value. The reason is that at the i-th iteration of the algorithm, 
the updated R^(t) is calculated on 

the basis of X®(t). However, although X®(t) = X (i_1) (t) + A (i_1) (t), it can happen 
that i2«(t) « R^(t) since ^{Fy^-^t) + A^t)]} « ^{F^X^t)]}. 
In this case, the subsequent iterations will be not able to further improve the result. 

Another problem was recently identified by Deodatis & Micaletti (2000): after the 
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first iteration the field X^(t) is no longer strictly Gaussian and therefore the mapping 
(27) will not give a R^(t) with the correct marginal distribution. These authors 
provide a modified version of the original algorithm of Yamazaki & Shinozuka where, 
after the first iteration, F x is substituted by an empirical distribution of jW(t). 
Actually, such a method works well also in case of very skewed distributions and/or 
non-smooth target power spectra. Unfortunately, it is very expensive with respect 
to computational time which makes its use problematic in practical situations (e.g. 
simulation of sizeable random fields) ; 

- in the case of homogeneous and isotropic iV-dimensional random fields, it is necessary 
to work in the iV-dimensional Fourier domain. Furthermore, the entire procedure 
must be restarted for each new simulation. 



In spite of these problems, the algorithm described in this section maintains a certain 
interest since, contrary to other techniques, it can be easily adapted for the simulation of 
vector- valued random fields (Popescu, Deodatis & Prevost 1998). 



4.4. Fixing the Mean and the Variance 

In all the methods presented in the previous sections, the mean and the variance of R(t) 
are fixed by F R (r). However, without modifying the correlation structure, one can force 
R(t) to have a given mean //* and a variance a\ by means of the following transformation 

fl.(t)=jx.+ i * (t) ~* i * a,. (44) 
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5. SOME POSSIBLE APPLICATIONS 

As already reported in the Introduction, the numerical simulation of non-Gaussian 
random fields can be used in understanding many both experimental and theoretical 
physical problems. 

Hot topics in Astrophysics and Cosmology, where the techniques described in this work 
find quick applications, can be easily identified as: 

• the simulation of continuous maps to match the properties of sky backgrounds. For 
instance, very deep maps of the extragalactic IR sky from space are plagued by the 
presence of Galactic Cirrus emission and at very small scales from source confusion; 

• the generation of non-Gaussian initials conditions for the iV-body simulations (see 
figure 8). Indeed, these conditions can be obtained by interpreting a non-Gaussian 
random fields as a density field. In reality, such approach is not new. However, in 
the past the initials conditions were simulated by using only specific distributions 
functions as, for example, the Lognormal (Moscardini et al. 1991; Coles & Jones 
1991) and the Chi-Square (Scoccimarro 2000) ones; 

• the reconstruction of the missing parts of experimental maps (e.g. angular distribution 
of i7L4 5" galaxies). Indeed, it is sufficient to transform the non-Gaussian field in 

a Gaussian one via the inverse of the mapping (27), to carry out the desired 
reconstruction through one among the many techniques available for the Gaussian 
case (e.g. Rybicki & Press (1992)), and then to trasform back the resulting field via 
the mapping (27). This techniques is described in detail in Sheth (1995) but, again, 
it is specialized to the Lognormal case. 
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6. A PRACTICAL EXAMPLE 

The approach presented in this paper is going to be used to simulate the FIR sky as 
it will be observed by the HERSCHEL Satellite (Pilbratt 2001). To mock extragalactic 
catalogues, built on the basis of theoretical modeling of the expected number of FIR sources, 
the emission from "local" (Galactic and interplanetary) backgrounds has to be added, in 
order to reproduce realistic observing conditions (Andreani et al. (2000) and Andreani 
et al. (2001)). Here, as an example, the reproduction of a typical Galactic background, 
starting from an observed sky region, is briefly outlined. 

Figure 9a shows a sky map, observed by the ISOPHOT camera (Lemke et al. 1996) on 
board of the ISO Satellite (Kessler et al. 1996) at 175/iin, with a projected size of roughly 
24' x 24' (Dole et al. 2001). The histogram of its values (see figure 9b) shows that the 
reproduction of such map requires the use of non-gaussian techniques. In particular, we 
have choosen the "change of variable method" with the following adjustments: 

- the PDF of the pixel values of the original image is estimated via the Johnson 
parametric method (for details see Vio et al. (1994)). With such an approach it is 
possible to build the mapping (27), as well as its inverse, in closed form. That results 
in a much less expensive numerical cost than in case of the use of the more popular 
histogram (see Figure 9b); 

- to infer the correlation function pxij), necessary for the numerical generation of 
the gaussian random field we prefer not to invert equation (6). Instead, we choose 
a set of forty-one equidistant values for px in the range [—0.2, 1.0], to determine 
the corresponding values of pn via the numerical integration of equation (6), and 
then interpolate the resulting points via a cubic-spline. In this way, again, a lot of 
computational effort is saved. The result is shown in figure 9c. 
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One of the possible simulations of the original map is shown in Figure 9d. 

7. SUMMARY AND CONCLUSIONS 

In this paper we have considered numerical simulation of non-Gaussian, scalar random 
fields R(t), with prescribed correlation structure Pr{t) and one-dimensional marginal 
probability distribution Fr, based on the transformation R(t) = g(X(t)) of a Gaussian 
random field X(t). In general, the definition of a function g(.), able to map the standard 
random Gaussian variable x in a random variable r with the required F R , is not a difficult 
task. Problems are found when the simulated fields has to have a desired correlation 
structure, since in general pxij) ^ Pr(t). The determination of the appropriate pxij) is 
achieved using various techniques. The most effective method is that providing a closed 
relationship between Px(t) and Pr(t). Unfortunately, this approach can be followed only 
in a very limited number of cases. Therefore, in the practical applications, very often it is 
necessary to resort to numerical techniques. 

Here, we have presented three approaches: the "change of variable method", the 
"Hermite expansion method" and the "method of Yamazaki & Shinozuka" . Whenever 
possible, the first one has to be adopted since, contrary to the other two, it is able to 
provide exact results. The only limitation concerning this method is that typically it 
requires the calculation of a large number of double integrals. In certain sitations that could 
be computationally too expensive. In this case, it is more convenient to use the "Hermite 
expansion method" since more robust and versatile than the "method of Yamazaki & 
Shinozuka" . This last method maintains a certain interest since it can be easily generalized 
for the numerical generation of vector-random fields. 

The authors warmly thank George Deodatis, Sabino Matarrese and Radu Popescu for 
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Fig. 1. — Examples of non-Gaussian random fields characterized by the same correlation 
function. 
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Fig. 2. — Relationship between the correlation function Pr(t) of some non-Gaussian random 
fields obtained via transformation (27) and the correlation function px(j) of the Gaussian 
fields used in such transformation. 
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Subroutine [a, Ng , px (t)] = HermCoeff[A r , fn, pn(r), e] 

1: u = RAND [N] 

2: x = F x 1 [u] 

3: r = F fl 1 [u] 

4 : f r .=0 

5 : a = MEAN[r] 

6 : r* = a 
7: k = 

8 : Do While : DIST[/ H , f r ,]>e 
9: k = k + l 

10: = HERMITE[/s,x] 

11 : a k = r T h k /(k\ N) 

12 : r* = r* + a k h k 

13 : / r , = EPDF[r*] 

14 : N H = k 

15 : End Do 

16 : Px (t) = IN VP OL [pn (r) , Ng] 

Return 



Subroutine R(t) = Field JHerm^jj, a, X(t)] 

1 : i?(t) = a 

2 : = 1 

3 : Do While : k < N H 

4 : R(t) = R(t) + a k ■ HERMITE[fc, X{t)\ 
5 : k = k+ 1 

6 : End Do 
Return 

RAND[.] = generator of uniform random numbers, MEAN[.] = mean value, 
DIST[., .] = distance measure, HERMITE[/s, .] = value of the k th Hermite poly- 
nomial given by equation (31) , EPDF[.] = empirical probability density function, 
INVPOL = polynomial inversion of equation (40). NB. Subroutine Field-Herm 
is a plain implementation of equation (32) with no claim of either numerical or 
computational efficiency. 
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Fig. 4. — Approximation of the \\ PDF obtained through the Hermite polynomial expansion 
with N H ranging from 1 to 4. For all the examples N = 16000. The quantity Gf is equal to 
\\fn — fr*\\ (see text), where f r * has been obtained via an histogram. 
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Subroutine R(t) = Field_IDF[0(k), S R (k), S , e] 

1: s£>(k) = 

2: sW(k) = S fl (k) 

3 : 5 = S 

4 : i = 1 

5 : Do While : DIST[5^(k), S(k)] > e 

6 : Jf'ljt) = IDFT [is'j'fk)! 1 / 2 e^< k) ] 

7: 1W = MEANpfM(t)] 

8: a xm =STD[X< i '(t)] 

9: l( i '(t) = [I (i) (t)-I (i) ]K,., 

10: i?«(t) =F^ 1 {F x [J?( < '(t)]} 

11: j?W = MEAN[i? (i) (t)] 

12: a fl(i) = STD[ J R< i >(t)] 

13: R (i) (t) = [JjW(t) - R (i) ]/a n(i) 

14: s£'(k) = |DFT[i?W(t)]| 2 

15 : If : DIST[s]| ) (k),5H(k)] < DIST^ -1 ' (k), S„(k)] 

16: k+ = FIND[k | s£'(k) = 0] 

17: S^(k+) = 1 

18: C(k) = S fl (k)/S«(k) 

19: k* =FIND[k | C(k) > 1 + d] 

20: C(k*) = l + <5 

21: 5<j +1) (k) = C(k)-5<j ) (k) 

22 : t = i + 1 

23 : Else : 

24 : S = SCALE[<5] 

25 : End If 

26 : End Do 

27 : R(t) = R (i - r > (t) + i?*^ 1 ' 

Return 



DIST[.,.]= distance measure, DFT= discrete Fourier transform, IDFT[.]= in- 
verse discrete Fourier transform, MEAN[.]= mean value, STD[.]= standard de- 
viation, FIND[b| condition]= it finds the elements of the array b satisfying 
"condition", SCALE[.]=scaling function. 
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Fig. 6. — Convergence rate of the Yamazaki & Shinozuka method concerning some type of 
non-Gaussian fields: Uniform (skewness= 0, kurtosis= —1.2), xi (skewness= 1.63, kurtosis= 
4) and Beta(4,2) (skewness= 0.47, kurtosis = 0.38). 
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Point processes obtained from the non-Gaussian random fields of figure 1. 
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Fig. 9. — a) Original FIR map obtained by the ISO satellite, b) Classical and Johnson 
histograms of the values of the map in the previous panel. The Gaussian PDF is plotted for 
reference, c) Correlation function of the original map and correlation function of the Gaussian 
map used in the numerical experiments (see text), d) Typical non-Gaussian simulated map. 



